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Abstract: 

A variational approach, based on a discrete representation of the chain, is used to calculate free 
energy and conformational properties in polyelectrolytes. The true bond and Coulomb potentials are 
approximated by a trial isotropic harmonic energy containing monomer-monomer force constants as 
variational parameters. By a judicious choice of representation and the use of incremental matrix 
inversion, an efficient and fast-convergent iterative algorithm is constructed, that optimizes the 
free energy. The computational demand scales as N 3 rather than N 4 as expected in a more naive 
approach. The method has the additional advantage that in contrast to Monte Carlo calculations 
the entropy is easily computed. An analysis of the high and low temperature limits is given. Also, 
the variational formulation is shown to respect the appropriate virial identities. The accuracy 
of the approximations introduced are tested against Monte Carlo simulations for problem sizes 
ranging from N = 20 to 1024. Very good accuracy is obtained for chains with unscreened Coulomb 
interactions. The addition of salt is described through a screened Coulomb interaction, for which 
the accuracy in a certain parameter range turns out to be inferior to the unscreened case. The 
reason is that the isotropic Gaussian Ansatz becomes less efficient with shorter range interactions. 

As a by-product a very efficient Monte Carlo algorithm was developed for comparisons, providing 
high statistics data for very large sizes - 2048 monomers. The Monte Carlo results are also used 
to examine scaling properties, based on low-T approximations to end-end and monomer-monomer 
separations. It is argued that the former increases faster than linearly with the number of bonds. 
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1 Introduction 



Polymers and polymer solutions play a profound role in our daily life, biologically and technolog- 
ically. This is of course one reason for the intense theoretical studies of polymers, but they have 
also because of their very nature been a challenge to theoreticians. Much theoretical work has been 
done in order to obtain a general understanding of neutral polymers, either in melts or in solution 
[1, 2]. Polyelectrolytes, on the other hand, have been less thoroughly investigated, despite their 
importance in many technical applications - for example in glue production or in food industry, for 
promoting flocculation or in pulp drying [3, 4, 5]. The term polyelectrolyte is sometimes used as a 
collective name for any highly charged aggregate. However, here we will restrict the meaning to a 
flexible molecule with several or many charged or chargeable sites; poly-L-glutamic acid, polyamines 
and polysaccharides are some typical representatives. 

The conformation of a flexible polyelectrolyte is a result of the competition between the covalent 
bonding forces, electrostatic interactions as well as more specific short ranged interactions. For 
example, poly-L-glutamic acid and several polysaccharides undergo helix to coil transition as a 
function of pH [6, 7]. This transition is obviously governed by electrostatic forces - similar structural 
transitions are seen for DNA [8]. Undoubtedly, both the polymer nature and the interaction between 
charged amino acids play an important role for the folding of a protein, as well as other solvent 
averaged forces. The present study should be seen as an attempt to approach the folding problem 
using what turns out to be a very powerful statistical mechanical variational technique. In the first 
step we will limit the study to linear polyelectrolytes in salt solution. 

Variational methods are standard techniques in quantum mechanics, but less so in statistical me- 
chanics, although variational principles were formulated many years ago [9]. One type of variational 
formulations starts off from an approximate free energy, which is optimized with respect to the par- 
ticle density [10]. For polymers a more fundamental approach is possible, by introducing variational 
parameters directly into the appropriate Hamiltonian. In the past this route has been followed in a 
number of polymer studies [11, 12, 13] and it has recently been revitalized by several groups [14, 15]. 
To the best of our knowledge, all these calculations have been concerned with continous chains and 
only one or at most a few variational parameters have been optimized. 

The present approach, of which some results were already published [16], is inspired by refs. 
[17, 18, 19]. It uses a discrete representation of the polymer, which not only allows us to investigate 
linear or cyclic polymers, but also polymers of arbitrary topology. Thus, e.g., a hyperbranched den- 
drite structure can easily be handled within this formalism. The number of variational parameters is 
proportional to N 2 , where N is the number of interacting units. In general, the variational approach 
is expected to be most accurate at high dimensions [17, 18]. Apparantly, this has discouraged the 
community from pushing the approach for three-dimensional polymers into a numerical confronta- 
tion. Also, using the method in a naive way would give scaling behaviour like N 4 , which would 
make the method less tractable for large sizes. We have also found empirically that such a naive 
implementation is plagued with bad convergence properties. In previous work [16] an algorithm 
was developed that lowers the computational costs to N 3 with controlled and nice convergence 
properties. 

The high and low T limits of the variational approach, which are accessible with analytical means, 
are derived in this paper. The low T expansions for various quantities yield results that very well 
approximate what emerges from the corresponding expansions in the exact theory. The fact that the 
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temperature of interest (room temperature) is fairly low on the temperature scale partly explains 
the success of initial numerical explorations of the variational approach [16] and motivates further 
studies. 

Besides being more realistic, the discrete chain also offers the possibility of direct comparison with 
Monte Carlo (MC) simulations, thereby giving an important indication of the accuracy. The final 
output of the present variational calculation is of course the minimized free energy, but one also 
obtains a matrix with all possible monomer-monomer correlations within the chain. This matrix is 
the starting point for the calculation of end-end and monomer-monomer separations and different 
kinds of angular correlations. 

MC simulations of flexible polyelectrolytes have only recently appeared in the literature and then 
limited to rather short chains [20, 21, 22, 23, 24, 25, 26, 27, 28], the longest chains being of the order 
of a few hundred monomers. So far simulations have dealt with both explicit representation of salt 
particles as well as implicit in the form of a screened Coulomb potential. The inclusion of short 
ranged interactions has also been studied as have the conformational changes associated with the 
titration of a flexible polyelectrolyte [27]. Most simulations have been carried out in the canonical 
ensemble, although the latter problem required the use of grand canonical MC simulations [26, 27]. 
The accuracy of the screened Coulomb potential has been investigated by several people and found 
to be an excellent approximation for not too high polyelectrolyte charge densities and in the absence 
of any multivalent ions [23, 29]. In order to obtain high statistics MC results a pivot algorithm and 
a recently developed hybrid scheme [30] were used. As a by-product these simulation results were 
also used to examine certain scaling properties, theoretically derived from an approximate low T 
expression. In contrast to the case of rigid bonds we find that the end-end separations scale faster 
than linearly with N. 

When confronting the variational approach with MC data the following results emerge in this work: 

• The variational approach has the unique property that it directly yields the free energy F. 
This is in contrast to MC simulations, where F is only indirectly accessible through elaborate 
integrations. For N=20, where comparisons are inexpensive, the variational free energy nicely 
agrees with MC results. 

• For unscreened Coulomb chains the success reported in ref. [16] survives to even larger systems 
(N=1024) for configurational properties like end-end correlations. Also for angular correlations 
and scaling behaviour the variational approach reproduces MC data very well. 

• When including salt through Debye screened Coulomb potentials the performance of the 
method deteriorates somewhat on the quantitative level. Even if the qualitative configura- 
tional picture agrees with that from MC data, the actual numbers for e.g. end-end correlations 
could differ up to 50 %. We attribute this discrepancy partly to the inability of Gaussians to 
reproduce short range interactions. 

One should also mention that other approximation schemes than the variational ones have been 
attempted in order to study polyelectrolyte conformations. In the mean field approximation it is 
only with the assumption of spherical symmetry that the mean field equations become tractable, 
but the results are not particularly encouraging [29]. In a cylindrical geometry, however, the mean 
field approximation behaves much more satisfactorily, but at the expense of large numerical efforts 
[31]. 
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This paper is organized as follows: In sect. 2 the basic Coulomb models (screened and unscreened) 
are presented. The variational approach is presented in sect. 3. The high and low T limits of 
the variational scheme are computed with analytical methods in sect. 4. A description of the MC 
method used in order to establish the quality of the variational method is found in Sect. 5. The 
results from confronting the variational approach with MC data are presented in sect. 6. Finally 
in sect. 7 a brief summary and outlook is given. Most of the detailed derivations are found in 
appendices - generics about the variational method (A), variational energies for a polyelectrolyte 
(B), virial identities (C), high and low T expansions (D) and zero temperature scaling properties 
(E). The disposition of the material into bulk text and appendices is such that the approach and 
results can be fully understood without reading the appendices. 



2 The Model 



2.1 The Unscreened Coulomb Chain 

In this model we consider a polyelectrolyte at infinite dilution and without any added salt. The 
polyelectrolyte counterions are thus neglected and the only electrostatic interactions are between 
the charged monomers. More explicitly, the polymer chain consists of N point charges connected 
by harmonic oscillator ("Gaussian") bonds. The potential energy for a chain then takes the form 



JV-l 



fei -«4&W^E^ (i) 

i=i i<j 1 J 1 

Here, the position of the ith charge, and 

x^j — x^ Xj (2) 

while q is the monomer charge and e r eo is the dielectric permittivity of the medium. We use the 
tilde notation E, iq, etc. for physical quantities in conventional units, and reserve E, x^, etc. for 
dimensionless ones, which will be used in the theoretical formalism below. 

The force constant can be reexpressed in terms of the iV = 2 equilibrium distance tq, given by 

„- ■>* 1 
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In terms of the dimensionless coordinates Xj, defined by 

X; = 7-oXi (4) 

the energy takes the form 

E = krt 



We will consider the system at a finite temperature T, which can be similarly rescaled, 

kr 2 
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with Jsb the Boltzmann constant. One then obtains for the Boltzmann exponent the simple expres- 
sion 

In other words, a polyelectrolyte at infinite dilution represents a two-parameter model where T and 
the number of monomers N are the only two non-trivial parameters of the system. 

Unless otherwise stated the following parameter values will be used throughout the paper; T=298K, 
e,.=78.3 and T"o=6A. Obviously, E depends only on the relative positions; the global center-of-mass 
position variable will have to be excluded from integrations over the coordinate space. 

The Gaussian and Coulomb energies are subject to a virial identity (see Appendix C), which in 
dimensionless units reads 

2(E G ) - (E c ) = 3{N - 1)T (8) 

Eq. (8) is a useful relation for checking the correctness and convergence behaviour in MC simulations 
and we find that it is in general obeyed to 0.3% or better. 



2.2 The Screened Coulomb Chain 



To treat a single polyelectrolyte in a solution at finite salt concentration becomes very costly in 
a MC simulation, since for reasonable salt concentrations the number of salt ions easily becomes 
prohibitively large, much larger than the number of polyelectrolyte monomers. The usual way to 
avoid this problem is to preaverage the degrees of freedom of the simple salt ions for some fixed 
configuration of the polyelectrolyte [32], thus defining salt averaged effective potentials - this is the 
basis of the electrostatic contribution in the classical DLVO potential [33]. In this way we may 
derive a screened Coulomb potential from a linearisation of the Poisson-Boltzmann equation for the 
salt. Eq. (1) is then replaced by 

where k is the Debye screening length for a 1:1 salt defined as 

/ 2N A c. 



In eq. (10) c s is the salt concentration in molars (M) and Na is the Avogadro's number. The 
Boltzmann factor will for the screened Coulomb potential contain the inverse dimensionless Debye 
screening length, re = r^k as an additional parameter, and with the parameter values given above 
we have for c s = 0.01 M, 0.1 M and 1.0 M the re-values 0.1992, 0.6300 and 1.992 respectively. Then 
eq. (7) is modified into 

K 1 / 1 ^ „ ^ P -K\Xij\ \ 
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The virial identity (see Appendix C) now takes the form 

2(E G ) - (E c ) - re(^V K l x «l) = 3{N - 1)T (12) 

2.3 Relative coordinates 

In the remainder of this paper, relative coordinates will be used; instead of the absolute monomer 
positions Xj, the bond vectors r^, 

ri = x i+1 -Xi, i=l,...,N -1 (13) 

will be used as the fundamental variables. In this way complications due to the translational zero- 
mode are avoided; in addition the convergence of the algorithm is considerably speeded up, especially 
at high temperatures. 

The energy of the screened Coulomb chain will then take the following form: 

N— 1 

1 x—T x e~ Kr " 
E(r) = E G + Ec = -J2* 2 i+J2 —T- ( 14 ) 

1 = 1 o 

where a runs over contiguous non-nil sub-chains, with 

*v = E r< ( 15 ) 

corresponding to the distance vector between the endpoints of the subchain. The unscreened chain 
results for re = 0. 



3 The Variational Approach 
3.1 The Gaussian Ansatz 

In refs. [18, 19, 16] the variational method of refs. [9, 17] (see Appendix A for a generic description) 
was revisited in the context of discrete chains of polyelectrolytes. The approach is based on an 
effective energy Ansatz Ey, given by 

Ev / T = \j2 G ~ii ^ - a *) • - ( 16 ) 

ij 

where defines an average bond vector, around which Gaussian fluctuations are given by the 
symmetric positive-definite correlation matrix Gij, the matrix inverse of which appears in the energy. 

Using this effective energy, the exact free energy F = —T log Z of the polymer is approximated from 
above [9] by the variational one 

F = F v + (E - E v ) v > F (17) 
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where Fy = —TlogZy, and ()v refers to averages with respect to the trial Boltzmann distribution 
exp(-E v /T). 



The parameters Gij and are to be determined such that the variational free energy F is minimized; 
we note that the number of variational parameters increases with N like N 2 . The resulting effective 
Boltzmann distribution is then used to approximate expectation values (/) by effective ones (f)v ■ 
Thus, we have e.g. 

(n) v = a; (18) 
(r; • Tj) v = Hi ■ a.j + 3Gij (19) 



For the polyelectrolyte systems treated in this study we will at high temperatures find a unique 
variational solution, characterized by a^ = 0; this defines a purely fluctuating solution. At low tem- 
peratures we find in addition a rigid solution with aligned a^ ^ 0. The latter is due to spontaneous 
symmetry-breaking; it ceases to exist at high temperatures, but will at low enough temperatures 
have the lower free energy. As discussed below, the rigid solution can typically be disregarded at 
normal temperatures. 

For potentials more singular than 1/r 2 , (E)v will be divergent, and the approach breaks down. 
However, such potentials are not physical and we do not consider this limitation of the approach a 
serious one. 

A non-trivial result of the scaling properties of the effective energy is that the virial identity, eqs. 
(8,12), will be respected by the above variational approach (see Appendix C). 



3.2 Using Local Fluctuation Amplitudes 

The minimization of F with respect to Gij and a^ gives rise to a set of matrix equations to be solved 
iteratively. These are considerably simplified, and the symmetry and positivity constraints on Gij 
are automatic, if Gij is expressed as the product of a matrix and its transpose: 

JV-l 
/i=l 



The interpretation of the local parameter Zi is simple - it is a fluctuation amplitude for the ith bond 
vector r^. We can write 

r; = ai + ^2 z if- J m ( 21 ) 

where each component of £ 1Z 3 is an independent Gaussian noise variable of unit variance. 
Similarly, we have for a subchain 

r (7 = ^r i Ea (7 + ^ z all J M , (22) 
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where a CT = X^go- a i anc ^ z " = Sign- z<- Thus, the noise amplitudes are additive. 
The matrix inverse of G can similarly be decomposed: 

G~i) = Wi • w,- (23) 
where Wi^ is the (transposed) matrix inverse of z^: 

Zj • Wj = % (24) 
Note that z^, and z CT are vectors, not in R 3 , but in R N_1 . 

The equations for a local extremum of F(a, z) are obtained by differentiation with respect to Zi and 

8F 8F 

^ = , — = 25 
"Zi oa^ 

3.3 The Unscreened Coulomb Chain 



In terms of and z^, the variational free energy for the pure Coulomb chain becomes, ignoring 
trivial additive constants (see Appendix B): 



26 



The equations for a minimum will be 



- = -STw.+Sz^-g-expj-^O 

(^2^) 



9F _ x - 

1 o-9i 



2 ^-erf| 

7T Z„ 



(27) 



(28) 



where the reciprocal vector is defined by eq. (24). 

These equations allow a purely fluctuating solution with a^ = 0. Setting a^ = the variational free 
energy simplifies to 



F = -3T log det z + \Y< Z * + Vf £ V a 



(29) 



which looks very much like the energy of an (N — l)-dimensional Coulomb chain with bonds z^, but 
with an extra entropy term (the first) preventing alignment of the ground state. The z derivatives 
become 

dF /2"v^ z CT , N 

— = -3T Wi + 3 Zi - J - V -f = 30 
dzi V 7T ^ 
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3.4 The Screened Coulomb Chain 



In the case of Debye screening the expression for F is modified to 
F = -3Tlogdetz + ^(3z 2 +a 2 ) 



2 

,2 



(31) 



where 

= exp(a; 2 /2) erfc(a;/v / 2) (32) 
The corresponding derivatives (cf. eq. (27)) take the form 



8F 
dzi 



8F 



-3Twi + 3zi (33) 
a; (34) 



Also here, a purely fluctuating solution is allowed. Setting = the variational free energy reduces 
to 



-3Tlogdet, + ^ + W^ 

i <r \ 



U (35) 



The z derivatives will be 



|^ = -3T Wi + 3 Zi - V 5| J - « 2 4) + K 3 Z 3 $ I = o (36) 

dZi f£ Zl \ V 7T j 



3.5 Implementation 

Due to the use of relative coordinates and of local noise amplitudes, a simple gradient descent 
method with a large step-size e can be used, that gives fast convergence to a solution of eqs. (25) 

8F 8F 
Azi = -e z —, Aa; = -e a — , (37) 
dzi den 

Further speed is gained by updating the reciprocal variables using incremental matrix inversion 
[16] - the increment in Wj due to Azi is given (exactly) by 

w,(w; • Az,) 

Aw,- = ^ (38) 

3 1 + Wi-Azi K 1 
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to be applied in parallel for j for fixed i. As a by-product the denominator (1 + • Az^) gives the 
multiplicative change in the determinant det z, needed to keep track of F. 

In order to maintain a reasonable numerical precision, the erf-related functions needed in the process 
are evaluated using carefully defined Taylor expansions or asymptotic expansions, depending on the 
size of the argument. 

As discussed above, the equations for a minimum are consistent with a purely fluctuating solution 
Hj = 0. Such a solution does indeed exist at all T; furthermore it is the only solution for high 
enough T. It turns out that for realistic choices of T one is in the region where this solution gives 
rise to good results (see figs. 1, 2 below). The additional solution with ^ appearing at low T is 
a symmetry-broken solution; to be realistic, such a solution should show an anisotropy also in the 
fluctuations, i.e. different amplitudes Zi for the fluctuations parallel and transverse to the direction 
defined by (the aligned) a^. This requires a more general Ansatz than eq. (16); theoretically, this 
will produce better low-T solutions, but for the Coulomb chain it leads to equations containing 
functions that are difficult to evaluate numerically, so we will not use this possibility in this paper. 
The incomplete symmetry-breaking partly explains the tendency for the a ^ solutions to produce 
inferior solutions. 

For the above reasons, and for reasons of continuity, we will in the numerical explorations use the 
a = solutions, where not otherwise stated. This also implies faster performance since only the 
z^-variables, with simplified updating equations, are needed. 

The complete algorithm will look as follows: 



1. Initialize Zi (and a^ if present) randomly, suitably in the neighborhood of a truncated high- 
or low-T series solution. 

2. For each i: 

• Update Zi (and a^) according to eqs. (37) with suitable step-sizes. 

• Correct all Wj according to eq. (38). 

3. Check if converged; if not, go to 2. 

4. Extract a^ and Gij = Zi ■ Zj, and compute variational averages of interest. 



Typical step-sizes are e z ps 1/6 and e a ps 1/2. The convergence check is done based on the rate 
of change and on the virial identity. The number of computations in each iteration step for this 
procedure is proportional to N 3 . The number of iterations required for convergence is a slowly 
growing function g(N). In total, thus, the execution time of the algorithm grows with N as N 3 g[N). 
In terms of CPU requirement convergence of a N = 40 chain (a^ = 0) requires 3 seconds on a DEC 
Alpha workstation. 
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4 High and Low T Results — Analytical Considerations 



Before embarking on a numerical evaluation of the variational approach with comparisons to MC 
results, it is interesting to see what can be gained from studying the high and low T limits, where 
analytical methods can be used. 

At the energy minimum, prevailing at T = 0, the polyelectrolyte will form a straight line. When the 
temperature is increased there will be a competition between the entropy and the repulsive Coulomb 
forces, and as T — > oo, the chain becomes Brownian, and the elongated or ordered structure is gone 
altogether. The temperature range of the transition from ordered to disordered structure is TV- 
dependent. As N increases at fixed T, the Coulomb force becomes relatively more important and 
the system effectively behaves as if T decreased: the polymer configuration becomes increasingly 
aligned. 

In the variational approach, as discussed above, the high temperature regime is characterized by a 
purely fluctuating solution, reflecting the Brownian nature of the chain at T — > oo. Such a solution 
survives as a local minimum also at lower T, where however also a rigid solution exists. Below a 
certain critical temperature T c , the latter gives the global minimum, indicating a first order phase 
transition. This is probably an artefact of the variational approach - in the MC simulations the 
system shows no evidence of possessing a phase transition (see section 5). The rigid solution mirrors 
the ordered elongated structure of the polymer at low T. 

With these qualitative arguments in mind, we turn to a more detailed investigation of the behaviour 
of the polyelectrolyte in the high and low T limits, together with an evaluation of the corresponding 
variational results. The unscreened and screened cases will be treated separately. 



4.1 The Unscreened Coulomb case 
4.1.1 High Temperature 

In the high T limit, the variational results can be expanded in 1/T (see Appendix D). Thus, for the 
expectation value of the Gaussian energy Eg, the first two terms of the expansion yield, 



which agrees with the exact result to the order shown; the first discrepancy occurs in the 0(T~ 2 ) 
term, as is in fact true for any quadratic expectation value (r^ry). By the virial identity, this also 
holds for the Coulomb energy (Ec) and for the total energy as well. 

For the individual bond lengths, the high T result can be written 




N-l 



(39) 




(40) 



where the last term is obtained from a continuum approximation, valid for large N . 
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4.1.2 Low Temperature 



In the low T limit (see Appendix D), the exact result for the total internal energy, expanded in 
powers of T, is, 

(E) = E + {3N - 5)T/2 + 0(T 2 ) (41) 

where Eq is the minimum energy at T = and 3N — 5 is the number of degrees of freedom, modified 
for the spherical symmetry (— > —2). It turns out that with the above first order low T correction, 
the energy, and thus also Eg, Eq and r mm (see below), are quite well approximated for the sizes 
and temperatures considered in this paper. 

At low T, the variational free energy is minimized by the rigid solution, for which the corresponding 
expansion is, 

(E) v =E + 3{N - l)T/2 + 0(T 2 ) (42) 

where the first term is the same as in eq. (41), while the second term is qualitatively correct for 
large N. 

Yet another low temperature expansion results from the purely fluctuating variational solution and 
it gives, 

(E) v = (6/tt) 1/3 Eo + 3{N - 2)T/2 + 0(T 2 ) (43) 

which shows a 24% discrepancy in the first term. The same factor, (6/V) 1 / 3 , results for any quadratic 
expectation value and for any single term in (Ec) in this limit. For r.m.s. distances, like the 
monomer- monomer distance r mm , 



i \" 2 



or the end-end distance r ee , 



this corresponds to an error of 11%. 



5>) ) ( 45 ) 



With the variational approach thus satisfactory in both temperature limits one can hope to find it 
a reasonable approximation also at finite temperatures, as will indeed be borne out in section 5. 



4.1.3 Zero Temperature 

Having low T expansions under control in terms of the T = configurational properties, the latter 
remain to be calculated. They are given by the minimum energy configuration, which is aligned, 

r; = bin (46) 

and unique (up to global translations and rotations). 
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The bond lengths bi > satisfy the equation 

* = £^ ( 4? ) 

where b a = ^2j ea .bj is the length of the subchain a. This equation cannot be solved analytically 
(except for very small N), but a fair large- N approximation can be obtained (see Appendix E for 
details). This gives for a distinct monomer-monomer bond the result 

1/3 

(48) 

As a consequence, we find that at zero temperature the average bond-length should scale logarith- 
mically with N , 

r mm oc {logNfl 3 (49) 

For the end-end separation this implies 

r ee oc N{\ogNY' 3 (50) 

This result is interesting, since it predicts a scaling faster than N; the extra logarithmic factor comes 
from the stretching of the harmonic bonds. 

The results above might seem as rather academic and of little practical impact. However, the 
low temperature expressions turn out to be quite accurate when compared with MC results. In 
other words ordinary room temperature and aqueous solution corresponds to a surprisingly "low" 
temperature. This indicates that the low temperature expansion might be a good starting point for 
further work in polyelectrolyte solutions. 



4.2 The Screened Coulomb case 

At high T, an analysis similar to the one carried out for the pure Coulomb chain, can be done 
for the screened chain (see Appendix D). Also there, the variational approximations to quadratic 
expectation values turn out correct to next-to-leading order; the same holds for E and Eq, while 
Ec (which is one order down) is correct to leading order. 

Also at low T, the results remain essentially the same for the screened chain. Thus, the rigid solution 
gives correct energies to lowest order in T, while the purely fluctuating solution does not. 



5 Monte Carlo Simulation Techniques 

For the numerical evaluation of the variational approach, the results were compared to those from 
MC simulations, which were performed in the canonical ensemble with the traditional Metropolis 
algorithm [34]. For short chains this is a straightforward procedure, but for linear chains consisting of 
more than about 100 monomers convergence problems appear. Typically for a chain of 40 monomers 
four million moves/monomer were required in order for the statistical fluctuations in the end-end 



\ « /T=0 



bi 



log const 

5 V N 
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separation to be less than one per cent. The energy terms and local conformational properties 
like monomer-monomer separations converged much faster. The addition of salt, i.e. the use of a 
screened potential, improves the convergence characteristics, but on the other hand its evaluation 
is more time-consuming than the pure inverse square root. Careful coding, with table look-ups for 
the inverse square root routine and, in particular, for the screened Coulomb potential, turned out 
to be more rewarding, reducing the computation time by almost a factor of three [35]. 

In order to treat longer chains with reasonable statistics we are forced to use more efficient algorithms 
like the pivot algorithm, first described in refs. [36, 37], with a high efficiency for linear chains on a 
lattice with short range interactions. Recently it has also been used successfully [38] for off-lattice 
simulations of a single polymer chain. It could be argued that the pivot algorithm should be even 
more efficient for chains with long range repulsive interactions like a charged polymer. The form 
of the pivot procedure used in this work can be described as a two step process, consisting of a 
random translation followed by a random rotation or vice versa. For a polymer chain with fixed 
bond lengths only random rotations will be used. The procedure is as follows: choose a monomer i 
and apply the same random translation to monomers i+1 to N . Then choose an axis at random and 
perform a random rotation of monomers i+1 to N around this axis. Evaluate the interaction energy 
between monomers 1 to i, and i+1 to N. This is a quadratic process in contrast to the single move 
algorithm, which only requires the evaluation of N pair interactions/move. Finally a Metropolis 
energy criterion is used to test for rejection or acceptance of the new configuration. We find that 
with a maximal random displacement of the order of 5-10 A and a maximal random rotation of tt, we 
reject approximately 50 % of the attempted moves. Typically, we generate 10 3 -10 4 passes (one pass 
= one attempted move/monomer) resulting in a statistical uncertainty in the end-end separation of 
approximately one per cent. The uncertainty in the average monomer- monomer separation and in 
the Coulomb and Gaussian energies is much less. Local averages, however, like the ith bond length, 
may have larger uncertainties, something that is also discussed by Madras and Sokal [37]. The 
pivot algorithm seems to be superior to the traditional single monomer procedure described initially 
already for chains with N > 20. We also have found that restricting the procedure to translational 
moves still makes it superior to the traditional algorithm. The pivot algorithm makes it feasible to 
simulate chains with more than one thousand interacting monomers. The major drawback with the 
pivot algorithm seems to be its limitations to linear chains or at least chains with simple topologies. 

Without excessive fine tuning of the translational and rotational displacement parameters, we find 
that the computational cost grows as N 3 . This power results from N 2 for each sweep of monomer 
moves, and an additional factor N from autocorrelations in quantities like end-end separations. 

Another efficient algorithm has recently been developed in ref. [30]. By identifying the slow modes 
in a Fourier analysis, one is able to use different random step lengths for different modes. This 
technique seems to be as efficient as the pivot algorithm and we have used it to check the accuracy 
of our simulations. For all cases investigated we obtain, within the statistical uncertainties, identical 
averages. The same is true for shorter chains where we also can use the original single monomer 
algorithm as a further test. 
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6 Numerical Results 



The superiority of the purely fluctuating variational solution over the rigid one, as discussed in 
section 5, is illustrated in figs. 1 and 2, where the variational results for r ee and r mm are compared 
to MC data for N = 20 and 80. 
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Figure 1: (a) r ee and (b) r mm as functions of T for an unscreened chain with N = 20. Filled 
circles represent MC data, and solid and dashed lines variational results, with a ^ and a = 0, 
respectively. 
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Figure 2: (a) r ee and (b) r mm as functions of T for an unscreened chain with N = 80. Same 
notation as in fig. 1. 



The remainder of this section will be devoted to comparisons of the variational approach to MC sim- 
ulation results focusing on (1) energies, in particular the free energy and (2) various configurational 
measures. The = variational solution will be consistently used. 
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6.1 Free and Internal Energies 



In table 1 the variational results for internal Coulombic and Gaussian energies are compared with the 
MC results; the relative deviations are seen to increase both with increasing c s and with increasing 
N . For fixed c s the deviations seem to converge to constant values at large N . Hence it should be 



N 


20 


40 


80 


160 


320 


512 


c s =0.0 E c (V) 
(MC) 


6.20 
5.25 


7.58 
6.30 


8.80 
7.28 


9.94 
8.16 


11.0 
8.99 


11.7 
9.53 


E G (V) 
(MC) 


6.65 
6.25 


7.40 
6.78 


8.08 
7.31 


8.66 
7.79 


9.20 
8.23 


9.54 
8.47 


c s =0.01 E c (V) 
(MC) 


3.55 
2.70 


3.80 
2.83 


3.95 
2.89 


4.02 
2.91 


4.05 
2.93 


4.07 
2.93 


E G (V) 
(MC) 


6.20 
5.70 


6.63 
6.00 


6.88 
6.15 


7.00 
6.22 


7.07 
6.26 


7.10 
6.27 


c.=0.1 E c (V) 
(MC) 


1.90 
1.35 


2.03 
1.38 


2.16 
1.40 


2.19 
1.41 


2.15 
1.42 


2.15 
1.42 


E G (V) 
(MC) 


5.40 
5.00 


5.68 
5.15 


5.86 
5.26 


5.94 
5.29 


5.93 
5.32 


5.94 
5.32 


c. = 1.0 E c (V) 
(MC) 


0.65 
0.40 


0.70 
0.40 


0.74 
0.41 


0.75 
0.43 


0.76 
0.42 


0.76 
0.42 


E G (V) 
(MC) 


4.35 
4.10 


4.53 
4.25 


4.61 
4.29 


4.67 
4.33 


4.69 
4.37 


4.70 
4.37 



Table 1: Average internal Coulombic and Gaussian energies per monomer in kJ/mol ■ monomer 
for unscreened and screened Coulomb potential. MC and V stands for Monte Carlo and variational 
calculations respectively. The salt concentration, c s , is given in molar (M). 

possible to extract "asymptotic" correction factors from comparisons at moderate N , do a variational 
calculation for a very large N , and predict what an MC calculation would give. 

A strong advantage of the variational approach is the direct access to the free energy, which is much 
more difficult to obtain in a MC simulation, requiring a cumbersome integration procedure. In 
order to evaluate the variational results we nevertheless attempt to estimate F(T) from MC data 
for N = 20 using the following procedure. In dimensionless units one has 

d{F/T) d , r _ EITj i 



;log 



J e - E /T dx = _l_ {E) (51) 



dT dT 

Thus, we can define an excess free energy with respect to some reference temperature T r as 

AF(T) = F(T)/T - F(T r )/T r = - f (E)dT/T 2 (52) 

which is then accessible in MC by a temperature integration of (E). In fig. 3 the excess free energy is 
shown as a function of T, for an N = 20 chain with T r corresponding to 1422 K. As can be seen the 
variational solutions for F(T) reproduce the extracted MC values very well in a wide temperature 
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Figure 3: The excess free energy AF(T) (see eq. (52)) from the MC data (filled circles) and from 
the variational approach (solid line with a ^ and dotted line with a = 0) for (a) c s = 0.0M and 
(b) c s = 0.1M respectively (N=10). 



interval. Comparisons for larger N are not feasible due to the indirect cumbersome MC extraction 
procedure. 

We remark that one possibly efficient alternative route to obtain free energies for different degrees 
of screening at fixed T in MC, would be to perform an integration in the Debye screening length and 
calculate the incremental excess free energy for a change Are. The advantage of such a procedure 
would be that every point along the integration path corresponds to a physically realistic situation, 
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N 


20 


40 


80 


160 


320 


512 


1024 


2048 




13.04 


13.60 


14.11 


14.57 


14.99 


15.26 


15.63 




MC 


12.56 


13.01 


13.43 


13.81 


14.17 


14.38 


14.68 


14.99 


diff. 


3.8% 


4.5% 


5.1% 


5.5% 


5.8% 


6.1% 


6.5% 




Tee V 


122 


277 


632 


1425 


3152 


5340 


11478 




MC 


119 


269 


606 


1347 


2958 


4985 


10651 


22507 


diff. 


2.5% 


3.9% 


4.3% 


5.8% 


6.6% 


7.1% 


7.8% 





Table 2: r mm and r ee in A for unscreened Coulomb potential (c s = 0.0M) as computed with the 
variational (V) and Monte Carlo (MC) methods. The errors originating from the MC runs are 
estimated to be O(0.2%). 

while the completely screened chain would serve as a reference state. 



6.2 The End-End Separation 

The end-end separation is a critical measure of the accuracy, being a global quantity with con- 
tributions from all bond-bond correlations. Unfortunately, it also turns out to be a complicated 
quantity from the convergence point of view in standard MC simulations, which means that it will 
have larger uncertainties than for example the average monomer-monomer separation - this is par- 
ticularly true for the pure Coulomb chain. Table 2 contains a comparison between variational and 
simulation results for r mm and r ee using the unscreened Coulomb potential. The result from the 
variational approach is impressive - the maximal deviation, 7.8% for N = 1024, from the MC results 
is well below the 11% bound discussed above. (Due to memory limitations on the local workstation 
N = 2048 has not been pursued with the variational approach.) This result can be compared to the 
results from a (spherical) mean field approach [39], which for the same system shows a deviation of 
20% already for N = 100. 

In section 4, r ee /N and r mm were conjectured to vary linearly with (log N) 1 / 3 for large N at zero 
temperature (cf. eqs. (49, 50)). In fig. 4, this linear dependence is tested on both MC and 
variational data at room temperature (298K), with a surprisingly good result. This indicates that 
room temperatures can be considered low for a reasonably long Coulomb chain with the chosen 
parameters. It is also clearly seen how the variational approximation to r ee exceeds the MC values; 
asymptotically we expect an 11% discrepancy, as discussed above. 

The conjectured zero-temperature scaling results seem to contradict the results by deGennes et al. 
[40] and Baumgartner [20] that r ee should scale linearly with N for an unscreened polyelectrolyte. 
However, the latter result is valid only if the monomer-monomer bonds are rigid, but with elastic 
bonds as in the present model, eq. (1), swelling is possible and r ee /N increases slowly with N due 
to the long-range Coulomb repulsion. 

Returning to the simulations by Baumgartner [20], we note that his effective temperature is almost 
a factor of ten larger than in the present study. Such a high temperature means that the chain 
behaviour is essentially brownian in character making it numerically difficult to detect the electro- 
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static expansion of the chain in a traditional MC simulation. The rigid monomer-monomer bonds 
used by Baumgartner also precludes the extra expansion predicted by eqs. (49,50). 
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Figure 4: (a) r mm as a function of (log N) 1 / 3 . Filled circles represent MC data and solid line the 
variational results. The dashed line is a linear fit to the MC data, (b) \ogr ee as a function of 
log(log-AT). Filled and open circles represent MC data and variational results respectively. The lines 
are linear fits. 



Table 3 contains the same quantities as table 2, but for a screened Coulomb potential. The agreement 
detoriates, when a small amount of salt is added. The screening reduces the Coulomb repulsion, 
which in a sense is the hard part in our variational calculation, and one would naively expect an 



19 



N 


20 


40 


80 


160 


320 


512 


C s — U.Ul 




V 


12.60 


12.87 


13.02 


13.10 


13.14 


13.16 






ML 


12.09 


12.24 


12.31 


12.34 


12.36 


12.37 






diff. 


4 2% 


5 1% 

\J . J. /U 


5 8% 

\J . (J /u 


6 1% 

\J . J. /U 


6 3% 

\J . tl /u 


6 4% 

U • i /U 




V 

1 ee 


V 


104 


201 


377 


680 


1188 


1710 






MC 


99.6 


183 


317 


521 


825 


1110 






diff. 


4.4% 


9.8% 


19% 


31% 


44% 


54% 


c s — U.l 




v 

V 


11.77 


11.90 


11.97 


12.01 


12.04 


12.04 






ML 


11.30 


11.35 


11.39 


11.38 


11.40 


11.40 






diff. 


4 2% 

±.* £j /U 


4 8% 

" . (J /(J 


5 1% 

\J . J. /U 


5 5% 

iJ . iJ /u 


5 6% 

iJ . \J /u 


5 6% 

iJ . \J /u 




Tee 


V 


78.2 


136 


231 


387 


640 


895 






MC 


72.9 


120 


192 


301 


459 


622 






diff. 


7.3% 


13% 


20% 


29% 


39% 


44% 


c, = 1.0 




V 


10.57 


10.69 


10.69 


10.69 


10.70 


10.70 






MC 


10.27 


10.29 


10.29 


10.30 


10.34 


10.32 






diff. 


2.9% 


3.9% 


3.9% 


3.8% 


3.5% 


3.7% 




^ee 


V 


55.0 


86.9 


137 


217 


343 


468 






MC 


52.1 


79.5 


122 


182 


283 


364 






diff. 


5.6% 


9.3% 


12% 


19% 


21% 


29% 



Table 3: r mm and r ee in A for the screened Coulomb potential (c s =0.01M,c s =0.1M and c s = 1.0M 
respectively) as computed with the variational (V) and Monte Carlo (MC) methods. The errors 
originating from the MC runs are estimated to be O(0.1%) 

improved accuracy with a decreased interaction. The opposite result is found and the discrepancy 
in the end-end separation becomes as large as about 40% with 10 mM of salt and N = 160. At 
sufficiently high salt concentration the agreement improves again, as it should, with the Coulomb 
repulsion completely screened. With a strong screening the Coulomb potential will have a short 
range, and for sufficiently large N the chain will be Brownian. This is also reflected in table 3, where 
the agreement for r ee is worst for intermediate chain lengths and improve again when N increases. 
The large discrepancy seen e.g. for c s =0.01 M and iV=160 might seem surprising, considering the 
excellent agreement found for the pure Coulomb chain, but it reflects the difficulty to properly 
emulate a short-range potential with a harmonic effective energy [18]. 



6.3 Monomer-Monomer Separations 

The variational results for the average monomer-monomer separation r mm is in excellent agreement 
with MC results for both the unscreened and screened cases - the largest error seen is of the order of 
5% (see tables 2 and 3). As for r ee the variational estimate is always larger than the MC value. This 
is also true for any single monomer-monomer separation (r?) 1 / 2 as can be seen in fig. 5. The shape 
of the curves are correctly reproduced by the variational solution and the largest discrepancy is not 
unexpectedly found in the middle of the chain, which may be explained by a stronger accumulated 
electrostatic repulsion there. These results can be compared to the mean field solution of ref. [25], 
in which gradually more and more pair interactions were treated explicitly and withdrawn from the 
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Figure 5: Bond lengths (r?) 1 / 2 along an N = 40 chain. Salt concentrations c s from top to bottom 
are M, 0.01M, 0.1M and 1 M respectively. Solid line represents variational results and dashed line 
Monte Carlo data. 



mean field. It was found that only after the inclusion of next-nearest neighbours, leading to lengthy 
numerical calculations, the curve shape became qualitatively correct. 

We discussed above the zero-temperature scaling for an unscreened chain, r mm oc (log .AT) 1 / 3 . When 
salt is added r mm increases much more slowly with chain length and for c s = l M, when r mm re -1 , 
it is essentially independent of chain length. One also notes that the individual bond lengths (rj) 1 / 2 
do not vary in the central part of the chain when salt is present (see fig. 5). 

Furthermore, for an unscreened chain, following eq. (48), we expect at T = the individual bond- 
lengths, raised to the power three, to vary linearly with u = log[s(l — s)], where s = i/N. In 
fig. 6a it is shown that this is approximately true, more so for the variational solution than for 
the MC results. The lower curve in fig. 6a, obtained with a screened Coulomb potential, shows a 
qualitatively different behaviour. Fig. 6b contains a similar graph for different N . 

The T = scaling relation for individual bonds, eq. (48), has the peculiar consequence that the 
length of a bond at the end of the chain becomes independent of N . This can be seen by rewriting 
the scaling relation as 



^/^{log^ + log^l-s)]} 



1/3 



[log i 



1/3 



(53) 



where the last expression holds for small i. Eq. (53) 
the first few bond lengths to be independent of N . 



also holds for the MC results, where we find 



21 



/ usr / users /bs / maims / common / fLg/polea.ps 



I usr / users /bs / manus / common / fig/pol6&.ps 



Figure 6: The cubed bond length (r|) 3 / 2 as a function of log(s(l — s)), with s = i/N the relative 
monomer position, (a) Upper line shows variational results and symbols MC data for N = 40, no 
screening. The zero temperature scaling corresponds to a straight line. The lower line shows MC 
data for a screened chain (c s = 0.01 M) for comparison, (b) Variational results for unscreened 
chains of varying size N . 

6.4 Angular Correlations 

In order to further test the variational solutions, we also have calculated angular correlations between 
bonds, 

C <,3 = -T== W 
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which roughly gives the average of the cosine between bonds. In fig. 7a, variational and MC data 
is shown for the neighbor correlation Cj^+i. It is seen that the variational Ansatz consistently 
overestimates the angle, more so in the presence of salt than without. Comparing figs. 5 and 7a, we 
find that the above discussed discrepancy between the MC and variational results for the end-end 
separation seems to be due to differences in angular correlations as well as in bond lengths. For the 
unscreened case the two sources seem to be of comparable magnitude, while for the screened case 
the angular correlations seem to be larger and the main cause of discrepancy. 
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Figure 7: Angular correlations (a) between neighbouring bonds, and (b) between the first and 
successive bonds. Solid lines represent variational and dashed lines MC results. The upper pair of 
curves are for an unscreened chain, while the lower pair is for c s = 0.01 M. (N = 80.) 
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Fig. 7b shows a more global angular correlation between the first and all successive bonds. The 
variational results for this quantity are in excellent agreement with MC data, both in the screened 
and in the unscreened chain. 



7 Summary and Outlook 



A deterministic variational scheme for discrete representations of polymer chains has been presented, 
where the true bond and Coulombic potentials are approximated with a trial isotropic harmonic 
energy. The variational parameters obey matrix equations, for which a very effective iterative 
solution scheme has been developed - the computational demand is N 3 . 

The high and low T properties of the variational approach has been analyzed with encouraging 
results. Also, the approach is shown to obey the relevant virial identities. 

In contrast to MC simulations, the free energy is directly accessible with the variational method. 

When confronting the results from the method with those from MC simulations, very good agreement 
is found for configurational quantities in the case of an unscreened Coulomb interaction (the error 
is within 11 %). 

In the screened case the method does not reproduce the MC results equally well although the 
qualitative picture of conformational properties is there. We attribute this problem to the difficulty 
for a Gaussian to emulate short range interactions. 

Recently, MC simulations were pursued for titrating Coulomb chains [41]. For such systems the 
Coulomb potential of eq. (1) is modified to 



where the binary variables Si are either 1 or depending whether monomer i is charged or not. 
Thus minimizing E now also includes a combinatorial problem - deciding where the charges should 
be located. Variational techniques related to the ones used in this paper have been successfully 
used in pure combinatorial optimization problems [42], where again tedious stochastic procedures 
are replaced by a set of deterministic equations. Along similar lines, the approach of this paper can 
be modified to allow for a variational treatment also of the titrating problem. 

The variational approach is also directly applicable to more general topologies - bifurcations pose 
no problems. Proteins could also be treated in this way provided the traditional Lennard-Jones 
potentials are replaced by forms that are less singular at the origin. 




(55) 
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Appendix A. The Variational Approach — Generalities 



In this appendix we discuss the generic variational approach that is used in this paper for the 
particular problem of a Coulomb chain. We consider a generic system, with dynamical variables x 
in some multi-dimensional state-space, and assume that a real energy function E(x) is given. 

For an arbitrary probability distribution P(x) in an arbitrary state space, the free energy with 
respect to an energy E(x) is generally defined as 

F = (E) — TS (Al) 

where S is the entropy, 

S = -(\ogP) (A2) 
and expectation values are defined with respect to P. Writing P(x) as 

P(x)=±eM-E v (x)/T) (A3) 

with 

Z = J dx exp(-E v {x)/T) (A4) 

the free energy can be written as 

F = -TlogZ + (E - E v ) (A5) 



Note that 



/\ / Ev — E 

exp(- J F/T) = Zexp/ " y 



< 



(A6) 

exp (^ Ev ~ E ^j ^ = J dx exp{-E/T) 



= exp{-F/T)\ Ev=E 

where the inequality is due to the convexity of the exponential function. Thus, F is bounded 
from below by its value for Ey = E, corresponding to the proper Boltzmann distribution, P(x) oc 
exp(-£(z)/T). 

The variation of F due to a variation bEy is given by 

8F = -T((8E V (E - E v )) - (8E V )(E - E v )) (A7) 
= -T(6E V {E - E v )) c 

where (ab)c stands for the connected expectation value (cumulant), (ab) — (a)(6). 

The idea of the variational approach is to choose a suitable simple Ansatz for the variational energy 
Ey, with a set of adjustable parameters en, i = 1, . . . , N p , the values of which are to be chosen so as 
to minimize the variational free energy. Demanding the vanishing of the variation of the free energy 
due to variations in the parameters en then leads to the general equations for an extremum: 

BE \ V 

-^(E-Ev)j =0, i=l,...,N p (A8) 
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where () v denotes an expectation value based on the variational Boltzmann distribution. This 
determines the optimal values of the parameters. Exact expectation values are then approximated 
by the corresponding variational ones. 
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Appendix B. The Variational Free Energy for the Chain 



In this appendix we derive the expressions for the variational free energy (eqs. (26, 31)), for the 
unscreened as well as the screened Coulomb chain, with or without translational parameters in the 
variational Ansatz. 



Unscreened Coulomb Chain 

For the specific case of the Coulomb chain of length N , the energy amounts to 



(Bl) 



fB2) 



where a is a contiguous subchain. 
Gaussian Parameters Only 

We consider first a pure Gaussian variational Boltzmann distribution, corresponding to 

E vi T = \Y. G ~i) v i- v i 

The parameter matrix G is forced to be symmetric and positive-definite by expressing it as 

JV-l 

Gij = z i ' z j = ^ ] z ifi z jfi (B3) 

The general expression for the variational free energy is 

F = -T\ogZ v - (E v )v + (E) v (B4) 

where the expectation values are with respect to the normalized variational Boltzmann distribution 
exp(—Ev/T)/Zv- The first two terms are trivial to compute. The first is 

3T 

-T\ogZ v = — logdetG -1 = -3Tlogdet.z (B5) 
apart from a trivial constant that can be neglected, as can the second term, 

-{E v )v = -\{N-l)T (B6) 



The last term, 



' a I V 
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consists in a sum of terms, each amounting to the variational expectation value of a simple function 
of a Gaussian vector variable r CT , of which is a special case. Its probability distribution is given 

by 

P(r CT ) oc exp (-^j (B8) 

with 

z<T = j2 z i ( B9 ) 

Thus, we have 

(*i)v = 3z? (BIO) 

and 

vX = ^ vfM (B11) 

Summing up, the variational free energy takes the form 

F = -3Tlogdet z + ^J2 z i + U iM ( B12 ) 

i (T 

to be minimized with respect to the variational parameters z^. 
Gaussian and Translational Parameters 

For the more general variational Ansatz with additional translational parameters a^, 

E v/T=^J2 G~l (r; - a;) ■ (r, - a,) 

i,3 



(B13) 



the main difference is a translation of the individual probability distributions of eq. (B8), which 
now read 

PMocexp (- (r %" z ^ )2 ) (B14) 

with 

a<T = j2 a i ( Bis ) 

This gives 

(p?)v = 3*? + a? (B16) 

and 

-\ = ± el f(^-)=U 2 c (z (T ,a (T ) (B17) 
r„ / v o,o VV2W 

The variational free energy becomes 

F = -3Tlogdet z + \ ^(3zf + a?) + £ U?{ Z<7 , o ff ) (B18) 
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Screened Coulomb Chain 

Next we consider the Debye-screened version of the Coulumb chain, with the energy 



Gaussian Parameters Only 

For the variational Ansatz with only Gaussian parameters, eq. (B2), we need eq. (BIO) and the 
expectation value 

The variational free energy then reads 

F = -3Tlogdet U ?M ( B21 ) 



Gaussian and Translational Parameters 



Finally, if for the screened Coulomb chain, eq. (B19), also translational parameters are used in Ey , 
eq. (B13), we will need the following result in addition to eq. (B16), 



/ exp(-KTy) 
\ r <r 



exp(/c 2 4/2) 



exp(— rea CT )erfc 



-J2z 



exp(rea CT )erfc ( — ^J- — - 



The variational free energy will read 

F = -3Tlogdet z + \ £)(3b? + a?) + U?{z a ,a a 



(B22) 



(B23) 
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Appendix C. The Virial Identity 



In this appendix we derive the virial identities and show that these are respected by the variational 
approach. 



Exact Virial Identity 

For any system described by a Boltzmann distribution 

P(x) = -|exp(-.E(x)/T) (CI) 

with x £ 1Z D , and E rising as a power for large |x|, we will have 

^ J V • (f (x) exp(-E/T))dx = (C2) 

for e.g. any polynomial f , due to the integrand being an exact divergence. This is equivalent to 

T(V • f ) = (f • WE) (C3) 
Thus, by varying f , we can obtain an infinite set of identities for the system. 

The virial identity results from the particular choice f = x; in its general form it reads 

(x • WE) = TD (C4) 
where D is the dimension of x-space (x • V is the scaling operator). 

This is particularly useful if the energy E is given by a sum of terms E a homogeneous in x, 

x • WE a = \ a E a (C5) 

in which case the virial identity takes the simple form 

J2*a(E a ) = TD (C6) 



This applies e.g. to the case of the unscreened polyelectrolyte. There the scaling operator is given, 
in relative coordinates, by • V r; , and we have Ac = 2 and Ac = —1; the virial identity thus 

reads 

2(E G ) - (E c ) = 3{N - 1)T (C7) 



Variational Virial Identity 

The virial identity is preserved by the variational approach under certain conditions, to be specified 
below. For the generic system above, minimizing the free energy 

F = F V + (E-E v )v (C8) 
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w.r.t. the parameters en of a variational energy Ey , leads to (cf. eq. (A8) in Appendix A) 

- ° <° 9 > 

Now, choosing f = x(£ l — Ey) in eq. (C3) with the variational Boltzmann distribution, we have 

DT(E - E v ) v + T(x • V(E - E v )) v = ((E - E v )x ■ VE V ) V (CIO) 
On the other hand, since the virial theorem holds for Ey, 

DT=(x-VE v ) v (Cll) 
Substituting this into eq. (CIO), we obtain 

T(x • V(E - E v )) v = ((E - E v )x ■ VE V )% (C12) 



If the set of parameters en of Ey is such (and this is the crucial condition), that the scaling operation 
on Ey can be written in terms of derivatives with respect to the parameters, i.e. if 

3E 

x-VSy = ^Gi(a)-^ (C13) 

i 1 

then the righthand side of eq. (C12) vanishes at the minimum, due to eq. (C9), and we are left 
with 

(x • VE) V = (x • VE V ) V = DT (C14) 

which is what we desired. 

Note that the derivation only relies on a local extremum of the free energy. Thus, for the polymer, 
the virial identity, eq. (C7), is respected by both the rigid (a ^ 0) and the purely fluctuating (a = 0) 
solutions. 
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Appendix D. High and Low T Expansions 



High T expansions 
Exact results 



At high T, the chain size will be large, and accordingly, the Gaussian term will dominate over the 
interaction V in the energy expression, 



with a summed over contiguous subchains. 



2- ' ' - '' )n 



It is then natural to attempt an expansion in the perturbation V. For an arbitrary expectation 
value, we have the perturbative expansion 

</> = </>° " ^(fV)° c + ^(fW)° c (D2) 

where ( )£. refers to connected expectation values (cumulants) in the unperturbed Boltzmann distri- 
bution. Due to the singular behaviour of v(r) for small r in the (screened or unscreened) Coulomb 
case, only the first few terms will be finite. This indicates that expectation values cannot be ex- 
panded in a pure power series in T, and that logarithmic corrections will occur after the first finite 
terms. 

We are interested in quadratic expectation values of the type (r^ • ry). These can be combined to 
give e.g. the rms end-to-end distance r ee , the gyration radius, and the Gaussian energy (Eg) (and 
thereby, in the pure Coulomb case, the interaction energy (Ec) by the virial identity). 

For the pure Coulomb chain, v(r) = 1/r, we get the perturbative expansion 

(Pi • p,-) = 3T% + ^ £ L- 3 ' 2 + 0(T~ 2 ) (D3) 

<r3i,j 

where a denotes a contiguous sub-chain containing the ith and the jth. bond, and L a its total 
number of bonds. This leads to 



E G ) = 3(N - l)T/2 + yj ^ £ L^ 2 + 0(T" 2 ) (D4) 



and 



r 2 ee = Z{N -l)T+^j—Y J L 1 J 2 + 0{T- 2 ) (D5) 



Similar results are obtained for the case of a screened Coulomb interaction, v(r) = e Kr /r, where 
the expansion of a quadratic expectation value gives 

(r, • p,-) = 3T% + J^^l E ^ 5/2 + °( T " 3 ) (° 6 ) 

<r3i,j 
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Variational results 



The corresponding variational results can also be expanded at high T (where a^ = 0), by expanding 
the variational solution around the unperturbed value, which can be chosen as \/T8i tl . The 
variational approximation to a quadratic expectation value, (rj • Vj)v = 3z^ • Zj, can then easily be 
expanded. The results thus obtained reproduce the exact results, eqs. (D3,D6), correctly to the 
order shown. 

We conclude that, independently of screening, the high T variational results are correct to next-to- 
leading order for E and Eq, and to leading order for Ec- 



Low T expansions 



Here we will treat only the pure Coulomb case in detail; most of the discussion applies also to the 
screened case. 



Exact results 



At low T, expectation values can be expanded around the configuration that minimizes the energy. 
This expansion is slightly complicated by the the rotational degeneracy of the minimum. The 
results can be expressed in terms of the classical configuration, which is given by a straight line 
configuration. 



Let bi be the the bond-lengths at the energy minimum. These have to be computed numerically, by 
solving the equation 

<T3i " 

where b a is the length of a subchain containing the ith bond. Note that the above equation can be 
written as a matrix equation: 

fc = £*i E ^r = £ 5 *A (D8) 

Thus, b is an eigenvector of the matrix B with a unit eigenvalue. Similarly, we can define a whole 
series of tensors: 



Ai 



J2^ = bi ( D1 °) 



B ^ = ( Dn ) 



y - 

<r3i,j,k " 



fD12) 
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They are all symmetric, and contracting either with bi gives the tensor of rank one less. 



In addition, we need two more matrices, related to 5, 

U = (1 + 25)" 1 (D13) 
V = P(1-B)~ 1 P (D14) 

where P denotes the projection matrix onto the subspace orthogonal to b, which is deleted by 1 — 5. 
In terms of these tensors, we have the quadratic expectation- values at low T: 

(n-Tj) = bibj+T f Uij + 2V ij + tjy* +3^2c k i m (biU jk + bjU ik ){U lm - V lm ) j +0(T 2 ) (D15) 

where the first two terms of the T coefficient are the naive contributions from the longitudinal and 
transverse fluctuations. The rest are corrections due to the rotational degeneracy of the T = 
configuration, which is also responsible for the transverse zero-modes (of 1 — 5). 

From this we obtain e.g the average Gaussian energy, 

(E G ) = 1/3E + T{3N/2 - 11/6) + 0(T 2 )q (D16) 
from which, using the virial identity, we obtain 

(E c ) = 2/3So - 2T/3 + 0(T 2 ) (D17) 

and 

(E) = E + T(3N -5)/2 + 0(T 2 ) (D18) 

where Eq = 3/2^ b\ is the exact energy at T = 0. In the last equation, the T-coefficient is, as it 
should, half the number of degrees of freedom, not counting the two rotational zero-modes (and the 
three translational ones already removed). 

In the screened case, similar results can be obtained. In particular, eq. (D18) remains valid, though 
with a different Eq. 



Variational results 



Similarly, the variational results can be expanded at low T. We have to distinguish between the 
two different solutions. 

For the purely fluctuating solution with a^ = 0, the variational free energy is, at T = 0, 

Fo = (E) v = 3 -^ + ] [lj2j- (D19) 

This is obviously just the energy of an (N — l)-dimensional version of the chain, with modified 
coefficients. It is minimized by the aligned configuration 




(D20) 
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where bi are given by eq. (D7), and n is some unit vector in N — 1 dimensions. For small but finite 
T, we have to add the entropy term, — 3Tlogdet{z}, to Fo. This forces the configuration out of 
alignment. The resulting configuration can be obtained as a low-T expansion around the T = 
solution. The first correction to z will be of order VT, and since the T = Zi are aligned, the 
matrix inverse diverges - it will go like 1 / \pT . 

For the total energy, the leading correction can be obtained as follows. The equation to solve is 

3T Wi = Vi F (b) (D21) 
For the T = solution, Vi-Fo(zo) = 0. Thus, to lowest order, 

3T Wi = Vi v j^o(z ) • dzj (D22) 

3 

The leading energy correction will be 

1 - 3T 

dF = -J2 ViVjFo^dzjdzi = — w i • dz i ( D23 ) 

ij i 

Now, dzi consists of an aligned part and a transverse part, both oc VT, while for the aligned 
part is oc 1, while the transverse part is oc l/\/T. Thus, the leading contribution to dpQ comes from 
the transverse part. Taking the trace of the tranverse part of the identity • Zj = 8ij, and noting 
that the transverse part of z sits entirely in dz, we have to leading order 

dF = 3T( ^" 2) = d(E) v (D24) 



Using the virial identity, which holds also for the variational expectation values, we get to first order 
in T: 



i)i/3# 0+ 3T (^ 2 ) +0 (T 2 ) (D25) 



(E)v = -„ . 

7T Z 

(E G )y = l(V 3 £o+ T(3j !,~ 4) +°( T2 ) (° 26 ) 

(E c )v = l{-) 1/3 E -T+O{T 2 ) (D27) 

O 7T 

where Eq is the exact T = energy. Note that already the zero-order results are off by a factor 
(fO 1 ^ 3 ~ 1-24, but that the correction to E is correct in the high N limit. 

Similar results are obtained for the screened Coulomb chain. Most of the general analysis leading 
to eq.(D24) still holds, and we have e.g. 

{E) v =E' + 3T( ^~ 2) + 0(T 2 ) (D28) 

with E' ^ E . 

For the symmetry-broken ^ solution, the T = configuration is instead given by 

a; = bin , Zi=0 (D29) 
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with hi as in eq. (D7), and n an arbitrary unit vector in 1Z 3 . The small T corrections are obtained 
from expressing F as 

F = -3Tlogdet{z} + (E{in + £)zi M J„))j (D30) 

/* 

where are uncorrelated standard Gaussian noise variables. Expanding in z, we obtain 

F = -3Tlogdet{z} + E(a) + ^ z< ' Z -?' V< ' V j E ( a ) + ■■■ ( D31 ) 

ij 

Because the Coulomb term satisfies Laplace' equation, this is just (provided a ^ 0) 

# = -3Tlogdet{z} + £'(a) + -J^ZiZi (D32) 

i 

and the variational free energy separates in z and a for small z. Thus, minimum is obtained for 

a^ = iibi (D33) 



and z satisfies 

which means 

The energies become 



3Twi + 3zi = (D34) 



ti ■ Zj = T6ij (D35) 



, ( 3T(N - 1) , 
(E) v = E +^- 2 (D36) 

(E G ) V = 1 -E 0+ 3T{N 2 - 1) (D37) 

{E c )v = \e (D38) 

which is correct to lowest order, and for large N qualitatively correct to first order (except for Ec)- 

Again, the screened chain lead to similar results; in particular, the T = energies will be the correct 
ones. 
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Appendix E. Zero Temperature Scaling Properties 



The T = configuration of a pure Coulombic chain cannot be obtained analytically, but must be 
computed numerically. However, an approximate calculation can be done. The equation for the 
bond lengths bi in the elongated ground state configuration is given by eq. (D7): 



* = £ij ( E1 ) 



By assuming that the bond length is locally approximately constant, this can be approximated by 

6 ^^EEV-^+1)- 2 (E2) 



b 2 

1 k = l l=i 



This can be rewritten as 

JV-l 



6? « (-1 + min(Z, N - i) + min(Z, i))/l 2 (E3) 



i=i 

This in turn can be approximated by an integral, leading to 

-,1/3 



bi 

Defining s = i/N, this amounts to 



, ( i( N 
log I const 



(E4) 



bi « (log(const Ns(l - s))) 1/3 (E5) 

Eqs. (E4,E5) give a quite accurate picture of the variation of the bond lengths at T = 0. They also 
imply, that the typical ground state bond length should swell roughly as (log N) 1 / 3 for large N. 
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